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Abstract 

The life of a tsunami is usually divided into three phases: the gen- 
eration (tsunami source), the propagation and the inundation. Each 
phase is complex and often described separately. A brief description 
of each phase is given. Model problems are identified. Their formula- 
tion is given. While some of these problems can be solved analytically, 
most require numerical techniques. The inundation phase is less doc- 
umented than the other phases. It is shown that methods based on 
Smoothed Particle Hydrodynamics (SPH) are particularly well-suited 
for the inundation phase. Directions for future research are outlined. 
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1 Introduction 

Given the broadness of the topic of tsunamis, our purpose here is to 
recall some of the basics of tsunami modeling and to emphasize some general 
aspects, which are sometimes overlooked. The life of a tsunami is usually 
divided into three phases: the generation (tsunami source), the propagation 
and the inundation. The third and most difficult phase of the dynamics of 
tsunami waves deals with their breaking as they approach the shore. This 
phase depends greatly on the bottom bathymetry and on the coastline type. 
The breaking can be progressive. Then the inundation process is relatively 
slow and can last for several minutes. Structural damages are mainly caused 
by inundation. The breaking can also be explosive and lead to the formation 
of a plunging jet. The impact on the coast is then very rapid. In very shallow 
water, the amplitude of tsunami waves grows to such an extent that typically 
an undulation appears on the long wave, which develops into a progressive 
bore Chanson [2005]. This turbulent front, similar to the wave that occurs 
when a dam breaks, can be quite high and travel onto the beach at great 
speed. Then the front and the turbulent current behind it move onto the 
shore, past buildings and vegetation until they are finally stopped by rising 
ground. The water level can rise rapidly, typically from to 3 meters in 90 
seconds. 

The trajectory of these currents and their velocity are quite unpredictible, 
especially in the final stages because they are sensitive to small changes in the 
topography, and to the stochastic patterns of the collapse of buildings, and to 
the accumulation of debris such as trees, cars, logs, furniture. The dynamics 
of this final stage of tsunami waves is somewhat similar to the dynamics of 
fiood waves caused by dam breaking, dyke breaking or overtopping of dykes 
(cf the recent tragedy of hurricane Katrina in August 2005). Hence research 
on fiooding events and measures to deal with them may be able to contribute 
to improved warning and damage reduction systems for tsunami waves in the 
areas of the world where these waves are likely to occur as shallow surge waves 
(cf the recent tragedy of the Indian Ocean tsunami in December 2004). 

Civil engineers who visited the damage area following the Boxing day 
tsunami came up with several basic conclusions. Buildings that had been 
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constructed to satisfy modern safety standards offered a satisfactory resis- 
tance, in particular those with reinforced concrete beams properly integrated 
in the frame structure. These were able to withstand pressure associated with 
the leading front of the order of 1 atmosphere (recall that an equivalent pres- 
sure p is obtained with a windspeed U of about 450 m/s, since p = 
By contrast brick buildings collapsed and were washed away. Highly porous 
or open structures survived. Buildings further away from the beach survived 
the front in some cases, but they were then destroyed by the erosion of the 
ground around the buildings by the water currents Hunt and Burgers [2005] . 

Section 2 provides a description of the tsunami source when the source 
is an earthquake. In Section 3, we review the equations that are often used 
for tsunami propagation. Section 4 provides a short discussion on the energy 
of tsunamis. Section 5 is devoted to the run-up and inundation of tsunamis. 
Finally directions for future research are outlined. 

2 Tsunami induced by near-shore earthquake 

The inversion of seismic data allows one to reconstruct the permanent 
deformations of the sea bottom following earthquakes. In spite of the com- 
plexity of the seismic source and of the internal structure of the Earth, sci- 
entists have been relatively successful in using simple models for the source 
Okada [1985]. A description of Okada's model follows. 

2.1 Introduction 

The fracture zones, along which the foci of earthquakes are to be found, 
have been described in various papers. For example, it has been suggested 
that Volterra's theory of dislocations might be the proper tool for a quan- 
titative description of these fracture zones Steketee [1958]. This suggestion 
was made for the following reason. If the mechanism involved in earthquakes 
and the fracture zones is indeed one of fracture, discontinuities in the dis- 
placement components across the fractured surface will exist. As dislocation 
theory may be described as that part of the theory of elasticity dealing with 
surfaces across which the displacement field is discontinuous, the suggestion 
seems reasonable. 

As commonly done in mathematical physics, it is necessary for simplicity's 
sake to make some assumptions. Here we neglect the curvature of the earth, 
its gravity, temperature, magnetism, non-homogeneity, and consider a semi- 
infinite medium, which is homogeneous and isotropic. We further assume 
that the laws of classical linear elasticity theory hold. 
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Several studies showed that the effect of earth curvature is negligible 
for shallow events at distances of less than 20° Ben-Mehanem et al. [1970], 
Ben-Menahem et al. [1969], McGinley [1969], Smylie and Mansinha [1971]. 
The sensitivity to earth topography, homogeneity, isotropy and half-space 
assumptions was studied and discussed recently Masterlark [2003]. The au- 
thor used a commercially available code, ABACUS, which is based on a finite 
element model (FEM). Six FEMs were constructed to test the sensitivity of 
deformation predictions to each assumption. The main conclusion is that the 
vertical layering of lateral inhomogeneity can sometimes cause considerable 
effects on the deformation fields. 

The usual boundary conditions for dealing with earth's problems require 
that the surface S of the elastic medium (the earth) shall be free from 
forces. The resulting mixed boundary-value problem was solved a century 
ago Volterra [1907]. Later, Steketee proposed an alternative method to solve 
this problem using Green's functions Steketee [1958]. 



2.2 Volterra's theory of dislocations 

In order to introduce the concept of dislocation and for simplicity's sake, 
this section is devoted to the case of an entire elastic space. The second 
reason is that in his original paper Volterra solved the problem in this case 
Volterra [1907]. 

Let O be the origin of a Cartesian coordinate system in an infinite elastic 
medium, Xi the Cartesian coordinates {i = 1,2,3), and ej a unit vector in 
the positive Xj-direction. A force F = Fe^ at O generates a displacement 
field u^{P, O) at point P, which is determined by the well-known Somigliana 
tensor ^ ^ 

Ui{P,0) = - — {6ikr^nn - arjk), with a = t— — 7-. (1) 
Svr/i A + 2/i 

In this relation 6ik is the Kronecker delta, A and /i are Lame's constants, 
and r is the distance from P to O. The coefficient a can be rewritten as 
a = 1/2(1 — u), where u is Poisson's ratio. Later we will also use Young's 
modulus E, which is defined as 

^^ /i(3A + 2/i) 
A + /i 

The notation r j means dr/dxi and the summation convention applies. 

The stresses due to the displacement field (1) are easily computed from 
Hooke's law: 



(2) 
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We find 

k ( r> /ON OlF I X^XjXfi; fl 5j^iXj -\- 6p;jXi ^ijXfi\ 

a,^[F,u) = -— + ^3 J ■ 

The components of the force per unit area on a surface element are denoted 
as follows: 

rrik k 

where the Vj-s, are the direction cosines of the normal to the surface element 
Sokolnikoff and Specht [1946]. A Volterra dislocation is defined as a surface 
S in the elastic medium across which there is a discontinuity Auj in the 
displacement fields of the type 

Awj = u\ — ul = Ui + flijXj, (3) 

fijj ^ji- (4) 

Equation (3) in which Ui and Qij are constants is the well-known Weingarten 
relation which states that the discontinuity Aui should be of the type of a 
rigid body displacement, thereby maintaining continuity of the components 
of stress and strain across S. 

The displacement field in an infinite elastic medium due to a dislocation 
of type (1) is then determined by Volterra's formula Volterra [1907] 

Ukiyi,y2,y3) ■=Uk{yi) = JJ ^^iT^dS. (5) 

s 

Once the surface S is given, the dislocation is essentially determined by 
the six constants Ui and Vlij. Therefore we also write 

Uk{yi) = j; jj 4(P,Q)z/,rf^+^ JJ{x,a^i{P,Q)-x,a^i{P,Q)}uidS, (6) 

where Qij takes only the values Q12, ^23, ^31- Following Volterra Volterra 
[1907] and Love Love [1944] we call each of the six integrals in (6) an ele- 
mentary dislocation. 

It is clear from (5) and (6) that the computation of the displacement field 
Uk{Q) is performed as follows: A force Fe^ is applied at Q, and the stresses 
o'ijiP, Q) that this force generates are computed at the points on E. In 

particular the components of the force on S are computed. After multipli- 
cation with prescribed weights of magnitude Auj these forces are integrated 
over S to give the displacement component in Q due to the dislocation on 
E. 
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2.3 Dislocations in elastic half-space 

When the case of an elastic half-space is considered, equation (5) remains 
valid, but we have to replace a^j by another tensor u^j. This can be explained 
by the fact that the elementary solutions for a half-space are different from 
Somigliana solution (1). 

The uOij can be obtained from the displacements corresponding to nuclei 
of strain in a half-space through relation (2). Steketee showed a method of 
obtaining the six ufj fields by using a Green's function and derived lJi2, which 
is relevant to a vertical strike-slip fault. Maruyama derived the remaining 
five functions Maruyama [1964]. 

It is interesting to mention here that historically these solutions were first 
derived in a straightforward manner by Mindlin Mindlin [1936], Mindlin and Cheng 
[1950], who gave explicit expressions of the displacement and stress fields for 
half-space nuclei of strain consisting of single forces with and without mo- 
ment. It is only necessary to write the single force results since the other 
forms can be obtained by taking appropriate derivatives. Their method con- 
sists in finding the displacement field in Westergaard's form of the Galerkin 
vector Westergaard [1935]. This vector is then determined by taking a linear 
combination of some biharmonic elementary solutions. The coefficients are 
chosen to satisfy boundary and equilibrium conditions. These solutions were 
also derived by Press in a slightly different manner Press [1965]. 



Figure 1: Coordinate system adopted in this study and geometry of the 
source model 

Here, we take the Cartesian coordinate system shown in Figure 1. The 
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elastic medium occupies the region z < and the x— axis is taken to be paral- 
lel to the strike direction of the fault. In this coordinate system, ul (xi, X2, x^, ^1,^2, ^3) 
is the ith component of the displacement at (xi, X2, X3) due to the jth direc- 
tion point force of magnitude F at (^1,^,^3). It can be expressed as follows 
Mindlin [1936], Okada [1985, 1992], Press [1965]: 

u](xi, X2, X3) = ul^{xi,X2, -X3)-ul^{xi, X2, X3)+ul^{xi, X2, X3)+X3ul(^{xi, X2, X3), 

(7) 

where 
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The first term in equation (7), ul^{xi, X2, X3), is the well-known Somigliana 
tensor, which represents the displacement field due to a single force placed 
at (^1,^2,^3) in an infinite medium Love [1944]. The second term also looks 
like a Somigliana tensor. This term corresponds to a contribution from an 
image source of the given point force placed at (^1,^2,-^3) in the infinite 
medium. The third term, m^^(xi, X2, X3), and ^^^(xi, X2, X3) in the fourth 
term are naturally depth dependent. When X3 is set equal to zero in equa- 
tion (7), the first and the second terms cancel each other, and the fourth term 
vanishes. The remaining term, m:^^(xi, X2, 0) reduces to the formula for the 
surface displacement field due to a point force in a half-space Okada [1985]: 
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In these formulas i?^ = (xi — + (x2 — ^2)^ + ■Cl- 
in order to obtain the displacements due to the dislocation we need to 
calculate the corresponding ^fc-derivatives of the point force solution (7) and 
to put it in Steketee-Volterra formula (5) 
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2.4 Finite rectangular source 

Now, let us consider a more practical problem. We define the elementary 
dislocations Ui, U2, and U3, corresponding to the strike-slip, dip-slip, and 
tensile components of an arbitrary dislocation. In Figure 1 each vector rep- 
resents the direction of the elementary faults. The vector D is the so-called 
Burger's vector, which shows how both sides of the fault are spread out: 
D = u+ - u". 

A general dislocation can be determined by three angles: the dip angle 
6 of the fault, the slip angle 6, and the angle (p between the fault plane and 
Burger's vector D. This situation is schematically described in Figure 2. 



O 



Free surface 



X 




Figure 2: Geometry of the source model and orientation of Burger's vector 
D 

For a finite rectangular fault with length L and width W occurring at 
depth d (Figure 2), the deformation field can be evaluated analytically by 
changing the variables and performing integration over the rectangle. This 
was done by several authors Chinnery [1963], Iwasaki and Sato [1979], Okada 
[1985, 1992], Sato and Matsu'ura [1974]. Here we give the results of their 
computations. The final results represented in compact form are listed below 
using Chinnery's notation || to represent the substitution 

/(e, ^)ll = f{x,p) - f{x,p- W) - fix -L,p) + fix -L,p-W). 
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Let us introduce the following notation: 

p = y cos 6 + dsinS, q = ysm6 — d cos S, 

y = r] cos 5 + g sin (5, (i = 77 sin 5 — g cos S, 

= ^^ + V^ + q^ = ^^ + f + d\ = ^2 ^ 

The quantities t/i, U2 and are linked to Burger's vector through the 
identities 

f/i = |D| cos(/)cos^, f/2 = |D| cos0sin^, U^ = \D\smcf). 
For a strike-slip dislocation, one has 

Ui f iq r . . 
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For a tensile fault dislocation, one has 



Ui 



U3 

U2 = — 



2% \R{R + 7]) 

—dq 
2tt \r(R+J) 



I3 sin^ 6 



sin S 



U3 



yq 



2ti \R{R + i) 



cos b 



_R{R + n) 

jq 

R{R + V) 



arctan 



— arctan 



qR 
qR 



Ii sin^ 6 



— h sin^ 5 



2.4 Finite rectangular source 



11 



parameter 


value 


LJl\J cllig,lc U 


1 o 


Fault depth d, km 


25 


Fault length L, km 


220 


Fault width W, km 


90 


Ui, m 


30 


Young modulus -E, GPa 


9.5 


Poisson's ratio u 


0.23 



Table 1: Parameter set used in Figures 3, 4, and 5. 
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Figures 3, 4, and 5 show the free-surface deformation after three elemen- 
tary dislocations. The values of the parameters are given in Table 1. 

The traditional approach for hydrodynamic modelers is indeed to use 
elastic models similar to the model just described with the seismic param- 
eters as input to evaluate details of the seafioor deformation. Then this 
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Figure 3: Dimensionless free-surface deformation z/a after dip-slip fault. 
Here a is |D| (30 m in the present application). 
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Figure 4: Dimensionless free-surface deformation z/a after strike-slip fault. 
Here a is |D| (30 m in the present application). 
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-600 -SOD 



Figure 5: Dimensionless free-surface deformation z/a after tensile fault. Here 
D = (0, 0, U-i) and a = U^. 
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deformation is translated to tlie initial condition of the evolution problem 
described in the next section. A few authors have solved the linearized wa- 
ter wave equations in the presence of a moving bottom Hammack [1973], 
Todorovska and Trifunac [2001]. 



3 Propagation of tsunamis 

The problem of tsunami propagation is a special case of the general water- 
wave problem. The study of water waves relies on several common assump- 
tions. Some are obvious while some others are questionable under certain 
circumstances. The water is assumed to be incompressible. Dissipation is 
not often included. However there are three main sources of dissipation for 
water waves: bottom friction, surface dissipation and body dissipation. For 
tsunamis, bottom friction is the most important one, especially in the later 
stages, and is sometimes included in the computations in an ad-hoc way. In 
most theoretical analyses, it is not included. 

A brief description of the common mathematical model used to study 
water waves follows. The horizontal coordinates are denoted by x and y, and 
the vertical coordinate by z. The horizontal gradient is denoted by 



dx dy ^ 

The horizontal velocity is denoted by 

u{x,y,z,t) = {u,v) 

and the vertical velocity by w{x,y,z,t). The three-dimensional flow of an 
inviscid and incompressible fluid is governed by the conservation of mass 

V.u+f^; = (8) 

oz 

and by the conservation of momentum 

Du Dw dp . , 

Pm=-^P^ Pm=-P'-^z 

In (9), p is the density of water (assumed to be constant throughout the fluid 
domain), g is the acceleration due to gravity and p{x,y,z,t) the pressure 
field. 

The assumption that the flow is irrotational is commonly made to analyze 
surface waves. Then there exists a scalar function (p{x, y, z, t) (the velocity 
potential) such that 

U = V0, w=—. 

oz 
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The continuity equation (8) becomes 

VV+|| = 0. (10) 

The equation of momentum conservation (9) can be integrated into Bernoulli's 
equation 

which is valid everywhere in the fluid. The constant Pq is a pressure of ref- 
erence, for example the atmospheric pressure. The effects of surface tension 
are not important for tsunami propagation. 



3.1 Classical formulation 

The surface wave problem consists in solving Laplace's equation (10) in a 
domain Q(t) bounded above by a moving free surface (the interface between 
air and water) and below by a fixed solid boundary (the bottom).^ The free 
surface is represented by F{x, y, z, t) := rj{x, y,t) — z = 0. The shape of the 
bottom is given hj z = —h{x,y). The main driving force is gravity. 

The free surface must be found as part of the solution. Two boundary 
conditions are required. The first one is the kinematic condition. It can be 
stated as DF/Dt = (the material derivative of F vanishes), which leads to 

r/t + V0- Vr/-0^ = at z = ri{x,y,t). (12) 

The second boundary condition is the dynamic condition which states that 
the normal stresses must be in balance at the free surface. The normal stress 
at the free surface is given by the difference in pressure. Bernoulli's equation 
(11) evaluated on the free surface z = r] gives 

<Pt + l\^<P\^ + l<Pl + 9V = at z = r]{x,y,t). (13) 

Finally, the boundary condition at the bottom is 

V</) ■ V/i + 0^ = atz = -h{x,y). (14) 

To summarize, the goal is to solve the set of equations (10), (12), (13) 
and (14) for ri{x,y,t) and (f){x,y,z,t). When the initial value problem is 

^Tlie surface wave problem can be easily extended to the case of a moving bottom. 
This extension may be needed to model tsunami generation if the bottom deformation is 
relatively slow. 
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integrated, tlie fields ri{x,y,0) and (j){x,y, z,0) must be specified at t = 0. 
The conservation of momentum equation (9) is not required in the solution 
procedure; it is used a posteriori to find the pressure p once rj and have 
been found. 

In the following subsections, we will consider various approximations of 
the full water-wave equations. One is the system of Boussinesq equations, 
that retains nonlinearity and dispersion up to a certain order. Another one 
is the system of nonlinear shallow-water equations that retains nonlinearity 
but no dispersion. The simplest one is the system of linear shallow-water 
equations. The concept of shallow water is based on the smallness of the ratio 
between water depth and wave length. In the case of tsunamis propagating 
on the surface of deep oceans, one can consider that shallow-water theory is 
appropriate because the water depth (typically several kilometers) is much 
smaller than the wave length (typically several hundred kilometers). 

3.2 Dimensionless formulation 

The derivation of shallow- water type equations is a classical topic. Two 
dimensionless numbers, which are supposed to be small, are introduced: 



where c? is a typical water depth, a a typical wave amplitude and i a typical 
wavelength. The assumptions on the smallness of these two numbers are 
satisfied for the Indian Ocean tsunami. Indeed the satellite altimetry obser- 
vations of the tsunami waves obtained by two satellites that passed over the 
Indian Ocean a couple of hours after the rupture occurred give an amplitude 
a of roughly 60 cm in the open ocean. The typical wavelength estimated 
from the width of the segments that experienced slip is between 160 and 240 
km Lay et al. [2005]. The water depth ranges from 4 km towards the west of 
the rupture to 1 km towards the east. These values give the following ranges 
for the two dimensionless numbers: 



1.5 X 10"^ < a < 6 X 10-^ 1.7 x 10~^ < /3 < 6.25 x 10"^ (16) 



The equations are more transparent when written in dimensionless variables. 
The new independent variables are 



where Cq = y/gd, the famous speed of propagation of tsunamis in the open 
ocean ranging from 356 km/h for a 1 km water depth to 712 km/h for a 4 




(15) 



X = ix, y = iy, z = dz, t = B/cq, 



(17) 
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km water depth. The new dependent variables are 

rj = afj, h = dh, cf) = galcf)/ Cq. (18) 

In dimensionless form, and after dropping the tildes, the equations become 

/?V'0 + 0.. = 0, (19) 

+ = aiz = -h{x,y), (20) 

/??7i + a/3V0 ■ V?7 = (pz at z = ari{x,y,t), (21) 

f3(f)t + la(3\V(j)\'^ + la(f)l + (37] = at z = ar]{x,y,t). (22) 



So far, no approximation has been made. In particular, we have not used 
the fact that the numbers a and j3 are small. 

3.3 Shallow-water equations 

When /3 is small, the water is considered to be shallow. The linearized 
theory of water waves is recovered by letting a go to zero. For the shallow 
water-wave theory, one assumes that (3 is small and expand in terms of (3: 

(j) = (f)o + (3(j)i + (3^(j)2 + ■■■ . 

This expansion is substituted into the governing equation and the boundary 
conditions. The lowest-order term in Laplace's equation is 

<Pozz = 0. (23) 

The boundary conditions imply that 0o = (po{x,y,t). Thus the vertical 
velocity component is zero and the horizontal velocity components are inde- 
pendent of the vertical coordinate z at lowest order. Let 0ox = u{x,y,t) and 
4>oy = v{x,y,t). Assume now for simplicity that the water depth is constant 
[h = 1). Solving Laplace's equation and taking into account the bottom 
kinematic condition yields the following expressions for 0i and (j)2'- 

Mx,y,z,t) = -^{1 + zfiu^ + Vy), (24) 
02(x,2/,z,t) = j-,il + z)^[{V\). + {V'v)y]. (25) 

The next step consists in retaining terms of requested order in the free- 
surface boundary conditions. Powers of a will appear when expanding in 
Taylor series the free-surface conditions around z = 0. For example, if one 
keeps terms of order a(3 and (3"^ in the dynamic boundary condition (22) and 
in the kinematic boundary condition (21), one obtains 

P<j)ot - l/3\ut, + vty) + /3r] + la(3{u^ + v^) = 0, (26) 
f3[r]t + a{u7], + vr]y) + {1 + a7]){u, + Vy)] = i/32[(V\), + (V^i;) J.27) 
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Differentiating (26) first with respect to x and then to respect to y gives a 
set of two equations: 

ut + a{uu^ + vv^) + r/x - \(i{utxx + vt^y) = 0, (28) 
vt + a{uuy + vvy) + T]y - l(3{utxy + vtyy) = 0. (29) 

The kinematic condition (27) can be rewritten as 

rjt + HI + ar])], + [v{l + ar])]y = |/?[(V\). + (V't;),]. (30) 

Equations (28)-(30) contain in fact various shallow-water models. The so- 
called fundamental shallow-water equations are obtained by neglecting the 
terms of order 

Ut + a{uUx + vUy) + rjx = 0, (31) 
vt + a{uvx + vvy) + r]y = 0, (32) 
r]t + [u{l + a7])]^ + [v{l + a7])]y = 0. (33) 

Recall that we assumed h to be constant for the derivation. Going back to an 
arbitrary water depth and to dimensional variables, the system of nonlinear 
shallow water equations reads 

Ut + uu^ + vuy + gr]^ = 0, (34) 
Vt + uvx + vvy + gr]y = 0, (35) 
r]t + Hh + 7])l + [v{h + r])]y = 0. (36) 

This system of equations has been used for example by Titov and Syno- 
lakis for the numerical computation of tidal wave run-up Titov and Synolakis 
[1998]. Note that this model does not include any bottom friction terms. To 
solve the problem of tsunami generation caused by bottom displacement, 
the motion of the seafioor obtained from seismological models Okada [1985] 
and described in Section 3 can be prescribed during a time to. Usually to is 
assumed to be small, so that the bottom displacement is considered as an in- 
stantaneous vertical displacement. This assumption may not be appropriate 
for slow events. 

The satellite altimetry observations of the Indian Ocean tsunami clearly 
show dispersive effects. The question of dispersive effects in tsunamis is 
open. Most propagation codes ignore dispersion. A few propagation codes 
that include dispersion have been developed Dalrymple et al. [2006]. A well- 
known code is FUNWAVE, developed at the University of Delaware over the 
past ten years Kirby et al. [1998]. Dispersive shallow water-wave models are 
presented next. 
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3.4 Boussinesq equations 

An additional dimensionless number, sometimes called the Stokes num- 
ber, is introduced: 

S = ^«l. (37) 

For the Indian Ocean tsunami, one finds 

0.24 < 5 < 46. (38) 

Therefore the additional assumption that S* ~ 1 may be realistic. 

In this subsection, we provide the guidelines to derive Boussinesq-type 
systems of equations Bona et al. [2002]. Of course, the variation of bathymetry 
is essential for the propagation of tsunamis, but for the derivation the water 
depth will be assumed to be constant. Some notation is introduced. The po- 
tential evaluated along the free surface is denoted by t) := (f){x, y, r], t). 
The derivatives of the velocity potential evaluated on the free surface are de- 
noted by := y, ?7, t), where the star stands for x, y, z or t. 
Consequently, (defined for * 7^ -z) and have different meanings. They 
are however related since 

= $(,) + <l>(^)r/* . 

The vertical velocity at the free surface is denoted by W{x^ y, t) := 4>z{x, y, r], t). 
The boundary conditions on the free surface (12) and (13) become 

r^t + V$ ■ Vr/ - W{1 + Vr] ■ Vr/) = 0, (39) 
$i + 5^+||V$p-|W^^(l + Vr/-Vr/) = 0. (40) 

These two nonlinear equations provide time-stepping for rj and $. In addi- 
tion, Laplace's equation as well as the kinematic condition on the bottom 
must be satisfied. In order to relate the free-surface variables with the bot- 
tom variables, one must solve Laplace's equation in the whole water column. 
In Boussinesq-type models, the velocity potential is represented as a formal 
expansion, 

00 

<P{x,y,z,t) = Y,<P^^\x,y,t)z'^. (41) 

n=0 

Here the expansion is about z = 0, which is the location of the free surface at 
rest. Demanding that </> formally satisfy Laplace's equation leads to a recur- 
rence relation between 0'^"^ and Let (po denote the velocity potential 
at 2; = 0, Uo the horizontal velocity at z = 0, and Wo the vertical velocity 
at 2; = 0. Note that c^o and Wo are nothing else than (/)*^°^ and (f)^^\ The 
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potential can be expressed in terms of 0o and Wo only. Finally, one obtains 
the velocity field in the whole water column {—h < z < t]) Madsen et al. 
[2003]: 

u{x,y,z,t) = cos{z'V)uo + sm{z'V)wo, (42) 
w{x,y,z,t) = cos{z'V)wo — sm{z'V)uo. (43) 

Here the cosine and sine operators are infinite Taylor series operators defined 

by 

oo 2n °° 2n+l 

n=0 ^ ' n=0 ^ ' 

Then one can substitute the representation (42)-(43) into the kinematic 
bottom condition and use successive approximations to obtain an explicit 
recursive expression for Wo in terms of Uq to infinite order in /iV. 

A wide variety of Boussinesq systems can been derived Madsen et al. 
[2003]. One can generalize the expansions to an arbitrary 2— level, instead of 
the z = level. The Taylor series for the cosine and sine operators can be 
truncated. Fade approximants can be used in operators a.t z = —h and/or 
at z = 0. 

The classical Boussinesq equations are more transparent when written in 
the dimensionless variables used in the previous subsection. We further as- 
sume that h is constant, drop the tildes, and write the equations for one spa- 
tial dimension (x). Ferforming the expansion about z = leads to the van- 
ishing of the odd terms in the velocity potential. Substituting the expression 
for into the free-surface boundary conditions evaluated at 2; = 1 + aT]{x, t) 
leads to two equations in rj and 0o with terms of various order in a and (3. 
The small parameters a and /3 are of the same order, while r] and 0^ as well 
as their partial derivatives are of order one. 



3.5 Classical Boussinesq equations 

The classical Boussinesq equations are obtained by keeping all terms that 
are at most linear in a or /3. In the derivation of the fundamental non- 
linear shallow-water equations (31)-(33), the terms in (3 were neglected. It 
is therefore implicitly assumed that the Stokes number is large. Since the 
cube of the water depth appears in the denominator of the Stokes number 
{S = a/ (3 = aP' it means that the Stokes number is 64 times larger in a 
1 km depth than in a 4 km depth! Based on these arguments, dispersion is 
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more important to tlie west of tlie rupture. Considering the Stokes number 
to be of order one leads to the following system in dimensional form^: 

Ut + UU^ + QT]^ - \h^Utxx = 0, (44) 

r]t + [u{h + 7])]^ - lh\^^^ = 0. (45) 

The classical Boussinesq equations are in fact slightly different. They are 
obtained by replacing u with the depth averaged velocity 

i-ri 

u dz. 

-h 



They read 



Ut + uu:^ + gr]^ - g/i utxx = 0, (46) 
Vt+Hh + r])l = 0. (47) 



A number of variants of the classical Boussinesq system were studied by 
Bona et al., who in particular showed that depending on the modeling of 
dispersion the linearization about the rest state may or may not be well- 
posed Bona et al. [2002]. 

3.6 Korteweg— de Vries equation 

The previous system allows the propagation of waves in both the positive 
and negative x— directions. Seeking solutions travelling in only one direction, 
for example the positive x— direction, leads to a single equation for r], the 
Korteweg-de Vries equation: 



rjt + Co (^1 + —j T]^ + -Cod T]^^^ = 0, (48) 

where d is the water depth. It admits solitary wave solutions travelling at 
speed V in the form 

r]{x, t)=a sech^ i^J^{x -Vt)^, with ^ = Cq (l + ^) . 

The solitary wave solutions of the Korteweg-de Vries equation are of elevation 
(a > 0) and travel faster than cq. Their speed increases with amplitude. Note 
that a natural length scale appears: 



4# 
3a 



^Equations (44) and (45) could have been obtained from equations (28) and (30). 
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For the Indian Ocean tsunami, it gives roughly (. = 377 km. It is of the order 
of magnitude of the wavelength estimated from the width of the segments 
that experienced slip. 

4 Energy of a tsunami 

The energy of the earthquake is measured via the strain energy released 
by the faulting. The part of the energy transmitted to the tsunami wave is 
less than one percent Lay et al. [2005]. They estimate the tsunami energy to 
be 4.2 X 10^^ J. They do not give details on how they obtained this estimate. 
However, a simple calculation based on considering the tsunami as a soliton 



The value for the integral is 4/3. The numerical estimate for E is close to 
that of Lay et al. (2005). Incidently, at this level of approximation, there 
is equipartition between kinetic and potential energy. It is also important 
to point out that a tsunami being a shallow water wave, the whole water 
column is moving as the wave propagates. For the parameter values used so 
far, the maximum horizontal current is 3 cm/s. However, as the water depth 
decreases, the current increases and becomes important when the depth be- 
comes less than 500 m. Additional properties of solitary waves can be found 
for example in Longuet-Higgins [1974]. 

5 Tsunami run-up 

The last phase of a tsunami is its run-up and inundation. Although in 
some cases it may be important to consider the coupling between fluid and 
structures, we restrict ourselves to the description of the fluid flow. The 
problem of waves climbing a beach is a classical one Carrier and Greenspan 
[1958]. The transformations used by Carrier and Greenspan are still used 
nowadays. The basis of their analysis is the one- dimensional counterpart of 
the system (34)-(36). In addition, they assume the depth to be of uniform 
slope: h = —a; tan 6*. Introduce the following dimensionless quantities, where 




gives for the energy 
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is a cliaracteristic lengtli^: 

X = ix, T] = if], U 



U, t 



After dropping tlie tildes, the dimensionless system of equations (34)-(36) 
becomes 

Ut + UUx + 7]^ = 0, 

rjt + [u{—x tan 9 + ri)]j; = 0. 

In terms of the variable c, these equations become 

Ut + uUx + 2cCx + tan 6* = 0, 
2ct + cux + 2uCx = 0. 



The equations written in characteristic form are 

{u + 2c + ttane) 



m + ^^^'^d-x 



d , , d 



iu — 2c + t tan( 



0, 
0. 



The characteristic curves C'^ and C as well as the Riemann invariants are 

dx 

: — = M + c, n + 2c + t tan 9 = r, 
at 



C- : ^ 
dt 



u 



ti — 2c + t tan^ = s. 



Next one can rewrite the hyperbolic equations in terms of the new variables 
A and a defined as follows: 



A 

2 
a 

4 



r + s) 



u -\-t tan 6, 



c. 



One obtains 



X a 



-(3r + s) - ttan6' 
4 



-(r + 3s) 
4 



t tan^ 



0. 



■^In fact there is no obvious characteristic length in this ideahzcd problem. Some authors 
simply say at this point that £ is specific to the problem under consideration. 
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The elimination of x results in the linear second-order equation for t 



Since u + ttan6' = A/2, u must also satisfy (49). Introducing the potential 
0((T, A) such that 



after integrating once. Two major simplifications have been obtained. The 
nonlinear set of equations have been reduced to a linear equation for m or 
and the free boundary is now the fixed line cr = in the (cr, A)— plane. The 
free boundary is the instantaneous shoreline c = 0, which moves as a wave 
climbs a beach. 

The above formulation has been used by several authors to study the 
run-up of various types of waves on sloping beaches Carrier et al. [2003], 
Tadepalli and Synolakis [1994], Tinti and Tonini [2005]. For example, it has 
been shown that leading depression A^-waves run-up higher than leading 
elevation iV-waves, suggesting that perhaps the solitary wave model may 
not be adequate for predicting an upper limit for the run-up of near-shore 
generated tsunamis. 

There is a rule of thumb that says that the run-up does not usually exceed 
twice the fault slip. Since run-ups of 30 meters were observed in Sumatra 
during the Boxing Day tsunami, the slip might have been of 15 meters or 
even more. 

Analytical models are useful, especially to perform parametric studies. 
However, the breaking of tsunami waves as well as the subsequent fioodings 
must be studied numerically. The most natural methods that can be used are 
the free surface capturing methods based on a finite volume discretisation, 
such as the Volume Of Fluid (VOF) or the Level Set methods, and the family 
of Smoothed Particle Hydrodynamics methods (SPH), applied to free-surface 
flow problems Gomez-Gesteira et al. [2005], Gomez-Gesteira and Dalrymple 
[2004], Monaghan [1994]. Such methods allow a study of flood wave dynam- 
ics, of wave breaking on the land behind beaches, and of the flow over rising 
ground with and without the presence of obstacles. This task is an essential 
part of tsunami modelling, since it allows the determination of the level of 
risk due to major flooding, the prediction of the resulting water levels in the 
flooded areas, the determination of security zones. It also provides some help 
in the conception and validation of protection systems in the most exposed 
areas. 




(49) 




one obtains the equation 



6 Direction for future research 



26 



6 Direction for future research 

A useful direction for future research in the dynamics of tsunami waves is 
the three-dimensional (3D) simulation of tsunami breaking along a coast. For 
this purpose, different validation steps are necessary. First more simulations 
of a two-dimensional (2D) tsunami interacting with a sloping beach ought 
to be performed. Then these simulations should be extended to the case of 
a 2D tsunami interacting with a sloping beach in the presence of obstacles. 
An important output of these computations will be the hydrodynamic load- 
ing on obstacles. The nonlinear inelastic behaviour of the obstacles may be 
accounted for using damage or plasticity models. The development of Boussi- 
nesq type models coupled with structure interactions is also a promising task. 
Finally there is a need for 3D numerical simulations of a tsunami interact- 
ing with a beach of complex bathymetry, with or without obstacles. These 
simulations will hopefully demonstrate the usefulness of numerical simula- 
tions for the definition of protecting devices or security zones. An important 
challenge in that respect is to make the numerical methods capable of han- 
dling interaction problems involving different scales: the fine scale needed for 
representing the damage of a flexible obstacle and a coarse scale needed to 
quantify the tsunami propagation. 
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